# ---------------------------
# R version: 4.3.2
# Code: 12/04/2023
# Title: 'Executive Action that Lasts'
# Authors: Kenneth Lowande, Michael Poznansky
# Summary: Replication code for producing all tables and figures. 
# Please report errors to: lowande@umich.edu
# ---------------------------


# install.packages(c('arm',doBy','ggplot2'))
require(arm)
require(doBy)
require(ggplot2)
# ---------------------------

load('.../orders.Rda')

##############################
# Model and Hypothesis Tests #
##############################

# model
model=glm(revoked_next_term~
            foreign+war+election+inflation+admchg+endterm+trend37+
            issuing_divgov*opposing_next*next_term_divgov*immutable, data=dt[dt$ceremonial==0,], family="binomial")

# set up simulation
set.seed(6658)
N=1000
s=sim(model,N)
betas=slot(s,'coef')

# variables from model
vars=names(model$coefficients)[2:length(model$coefficients)]

######
# H1 #
######

# negative case
zero=dt[complete.cases(dt[,vars[1:11]]),vars[1:11]]
zero[,'immutable']=0

zero[,'issuing_divgov:opposing_next']=zero$issuing_divgov*zero$opposing_next
zero[,'issuing_divgov:next_term_divgov']=zero$issuing_divgov*zero$next_term_divgov
zero[,'opposing_next:next_term_divgov']=zero$opposing_next*zero$next_term_divgov
zero[,'issuing_divgov:immutable']=zero$issuing_divgov*zero$immutable
zero[,'opposing_next:immutable']=zero$opposing_next*zero$immutable  
zero[,'next_term_divgov:immutable']=zero$next_term_divgov*zero$immutable  
zero[,'issuing_divgov:opposing_next:next_term_divgov']=zero$issuing_divgov*zero$opposing_next*zero$next_term_divgov
zero[,'issuing_divgov:opposing_next:immutable']=zero$issuing_divgov*zero$opposing_next*zero$immutable  
zero[,'issuing_divgov:next_term_divgov:immutable']=zero$issuing_divgov*zero$next_term_divgov*zero$immutable  
zero[,'opposing_next:next_term_divgov:immutable']=zero$opposing_next*zero$next_term_divgov*zero$immutable  
zero[,'issuing_divgov:opposing_next:next_term_divgov:immutable']=zero$issuing_divgov*zero$opposing_next*zero$next_term_divgov*zero$immutable

# positive case
one=dt[complete.cases(dt[,vars[1:11]]),vars[1:11]]
one[,'immutable']=1

one[,'issuing_divgov:opposing_next']=one$issuing_divgov*one$opposing_next
one[,'issuing_divgov:next_term_divgov']=one$issuing_divgov*one$next_term_divgov
one[,'opposing_next:next_term_divgov']=one$opposing_next*one$next_term_divgov
one[,'issuing_divgov:immutable']=one$issuing_divgov*one$immutable
one[,'opposing_next:immutable']=one$opposing_next*one$immutable  
one[,'next_term_divgov:immutable']=one$next_term_divgov*one$immutable  
one[,'issuing_divgov:opposing_next:next_term_divgov']=one$issuing_divgov*one$opposing_next*one$next_term_divgov
one[,'issuing_divgov:opposing_next:immutable']=one$issuing_divgov*one$opposing_next*one$immutable  
one[,'issuing_divgov:next_term_divgov:immutable']=one$issuing_divgov*one$next_term_divgov*one$immutable  
one[,'opposing_next:next_term_divgov:immutable']=one$opposing_next*one$next_term_divgov*one$immutable  
one[,'issuing_divgov:opposing_next:next_term_divgov:immutable']=one$issuing_divgov*one$opposing_next*one$next_term_divgov*one$immutable

# simulate predicted probabilities
zero=cbind(1,zero)
one=cbind(1,one)
p.zero=apply((1/(1+exp(-as.matrix(zero)%*% t(betas)))),1,as.vector) 
p_zero=apply(p.zero,1,mean)
p.one=apply((1/(1+exp(-as.matrix(one)%*% t(betas)))),1,as.vector) 
p_one=apply(p.one,1,mean)
diff=p_one-p_zero

# results
h1_lb <- quantile(diff,.025)
h1_est <- mean(diff)
h1_ub <- quantile(diff,.975)

######
# H2 #
######

# effect of change in party of president for mutable orders

# negative case
zero=dt[complete.cases(dt[,vars[1:11]]),vars[1:11]]
zero[,'immutable']=0
zero[,'opposing_next']=0

zero[,'issuing_divgov:opposing_next']=zero$issuing_divgov*zero$opposing_next
zero[,'issuing_divgov:next_term_divgov']=zero$issuing_divgov*zero$next_term_divgov
zero[,'opposing_next:next_term_divgov']=zero$opposing_next*zero$next_term_divgov
zero[,'issuing_divgov:immutable']=zero$issuing_divgov*zero$immutable
zero[,'opposing_next:immutable']=zero$opposing_next*zero$immutable  
zero[,'next_term_divgov:immutable']=zero$next_term_divgov*zero$immutable  
zero[,'issuing_divgov:opposing_next:next_term_divgov']=zero$issuing_divgov*zero$opposing_next*zero$next_term_divgov
zero[,'issuing_divgov:opposing_next:immutable']=zero$issuing_divgov*zero$opposing_next*zero$immutable  
zero[,'issuing_divgov:next_term_divgov:immutable']=zero$issuing_divgov*zero$next_term_divgov*zero$immutable  
zero[,'opposing_next:next_term_divgov:immutable']=zero$opposing_next*zero$next_term_divgov*zero$immutable  
zero[,'issuing_divgov:opposing_next:next_term_divgov:immutable']=zero$issuing_divgov*zero$opposing_next*zero$next_term_divgov*zero$immutable

# positive case
one=dt[complete.cases(dt[,vars[1:11]]),vars[1:11]]
one[,'immutable']=0
one[,'opposing_next']=1

one[,'issuing_divgov:opposing_next']=one$issuing_divgov*one$opposing_next
one[,'issuing_divgov:next_term_divgov']=one$issuing_divgov*one$next_term_divgov
one[,'opposing_next:next_term_divgov']=one$opposing_next*one$next_term_divgov
one[,'issuing_divgov:immutable']=one$issuing_divgov*one$immutable
one[,'opposing_next:immutable']=one$opposing_next*one$immutable  
one[,'next_term_divgov:immutable']=one$next_term_divgov*one$immutable  
one[,'issuing_divgov:opposing_next:next_term_divgov']=one$issuing_divgov*one$opposing_next*one$next_term_divgov
one[,'issuing_divgov:opposing_next:immutable']=one$issuing_divgov*one$opposing_next*one$immutable  
one[,'issuing_divgov:next_term_divgov:immutable']=one$issuing_divgov*one$next_term_divgov*one$immutable  
one[,'opposing_next:next_term_divgov:immutable']=one$opposing_next*one$next_term_divgov*one$immutable  
one[,'issuing_divgov:opposing_next:next_term_divgov:immutable']=one$issuing_divgov*one$opposing_next*one$next_term_divgov*one$immutable

# simulate predicted probabilities
zero=cbind(1,zero)
one=cbind(1,one)
p.zero=apply((1/(1+exp(-as.matrix(zero)%*% t(betas)))),1,as.vector) 
p_zero=apply(p.zero,1,mean)
p.one=apply((1/(1+exp(-as.matrix(one)%*% t(betas)))),1,as.vector) 
p_one=apply(p.one,1,mean)
diff=p_one-p_zero

diff_mutable_h2 <- diff

# results
h2_lb_a <- quantile(diff,.025)
h2_est_a <- mean(diff)
h2_ub_a <- quantile(diff,.975)

# effect of change in party of president for immutable orders

# negative case
zero=dt[complete.cases(dt[,vars[1:11]]),vars[1:11]]
zero[,'immutable']=1
zero[,'opposing_next']=0

zero[,'issuing_divgov:opposing_next']=zero$issuing_divgov*zero$opposing_next
zero[,'issuing_divgov:next_term_divgov']=zero$issuing_divgov*zero$next_term_divgov
zero[,'opposing_next:next_term_divgov']=zero$opposing_next*zero$next_term_divgov
zero[,'issuing_divgov:immutable']=zero$issuing_divgov*zero$immutable
zero[,'opposing_next:immutable']=zero$opposing_next*zero$immutable  
zero[,'next_term_divgov:immutable']=zero$next_term_divgov*zero$immutable  
zero[,'issuing_divgov:opposing_next:next_term_divgov']=zero$issuing_divgov*zero$opposing_next*zero$next_term_divgov
zero[,'issuing_divgov:opposing_next:immutable']=zero$issuing_divgov*zero$opposing_next*zero$immutable  
zero[,'issuing_divgov:next_term_divgov:immutable']=zero$issuing_divgov*zero$next_term_divgov*zero$immutable  
zero[,'opposing_next:next_term_divgov:immutable']=zero$opposing_next*zero$next_term_divgov*zero$immutable  
zero[,'issuing_divgov:opposing_next:next_term_divgov:immutable']=zero$issuing_divgov*zero$opposing_next*zero$next_term_divgov*zero$immutable

# positive case
one=dt[complete.cases(dt[,vars[1:11]]),vars[1:11]]
one[,'immutable']=1
one[,'opposing_next']=1

one[,'issuing_divgov:opposing_next']=one$issuing_divgov*one$opposing_next
one[,'issuing_divgov:next_term_divgov']=one$issuing_divgov*one$next_term_divgov
one[,'opposing_next:next_term_divgov']=one$opposing_next*one$next_term_divgov
one[,'issuing_divgov:immutable']=one$issuing_divgov*one$immutable
one[,'opposing_next:immutable']=one$opposing_next*one$immutable  
one[,'next_term_divgov:immutable']=one$next_term_divgov*one$immutable  
one[,'issuing_divgov:opposing_next:next_term_divgov']=one$issuing_divgov*one$opposing_next*one$next_term_divgov
one[,'issuing_divgov:opposing_next:immutable']=one$issuing_divgov*one$opposing_next*one$immutable  
one[,'issuing_divgov:next_term_divgov:immutable']=one$issuing_divgov*one$next_term_divgov*one$immutable  
one[,'opposing_next:next_term_divgov:immutable']=one$opposing_next*one$next_term_divgov*one$immutable  
one[,'issuing_divgov:opposing_next:next_term_divgov:immutable']=one$issuing_divgov*one$opposing_next*one$next_term_divgov*one$immutable

# simulate predicted probabilities
zero=cbind(1,zero)
one=cbind(1,one)
p.zero=apply((1/(1+exp(-as.matrix(zero)%*% t(betas)))),1,as.vector) 
p_zero=apply(p.zero,1,mean)
p.one=apply((1/(1+exp(-as.matrix(one)%*% t(betas)))),1,as.vector) 
p_one=apply(p.one,1,mean)
diff=p_one-p_zero

diff_immutable_h2 <- diff

# results
h2_lb_b <- quantile(diff,.025)
h2_est_b <- mean(diff)
h2_ub_b <-quantile(diff,.975)

h2 <- diff_mutable_h2-diff_immutable_h2

# discussed in-text
quantile(h2,.025)
quantile(h2,.05)
mean(h2)
quantile(h2,.95)
quantile(h2,.975)

######
# H3 #
######

# negative case
zero=dt[complete.cases(dt[,vars[1:11]]),vars[1:11]]
zero[,'immutable']=0
zero[,'issuing_divgov']=0

zero[,'issuing_divgov:opposing_next']=zero$issuing_divgov*zero$opposing_next
zero[,'issuing_divgov:next_term_divgov']=zero$issuing_divgov*zero$next_term_divgov
zero[,'opposing_next:next_term_divgov']=zero$opposing_next*zero$next_term_divgov
zero[,'issuing_divgov:immutable']=zero$issuing_divgov*zero$immutable
zero[,'opposing_next:immutable']=zero$opposing_next*zero$immutable  
zero[,'next_term_divgov:immutable']=zero$next_term_divgov*zero$immutable  
zero[,'issuing_divgov:opposing_next:next_term_divgov']=zero$issuing_divgov*zero$opposing_next*zero$next_term_divgov
zero[,'issuing_divgov:opposing_next:immutable']=zero$issuing_divgov*zero$opposing_next*zero$immutable  
zero[,'issuing_divgov:next_term_divgov:immutable']=zero$issuing_divgov*zero$next_term_divgov*zero$immutable  
zero[,'opposing_next:next_term_divgov:immutable']=zero$opposing_next*zero$next_term_divgov*zero$immutable  
zero[,'issuing_divgov:opposing_next:next_term_divgov:immutable']=zero$issuing_divgov*zero$opposing_next*zero$next_term_divgov*zero$immutable

# positive case
one=dt[complete.cases(dt[,vars[1:11]]),vars[1:11]]
one[,'immutable']=0
one[,'issuing_divgov']=1

one[,'issuing_divgov:opposing_next']=one$issuing_divgov*one$opposing_next
one[,'issuing_divgov:next_term_divgov']=one$issuing_divgov*one$next_term_divgov
one[,'opposing_next:next_term_divgov']=one$opposing_next*one$next_term_divgov
one[,'issuing_divgov:immutable']=one$issuing_divgov*one$immutable
one[,'opposing_next:immutable']=one$opposing_next*one$immutable  
one[,'next_term_divgov:immutable']=one$next_term_divgov*one$immutable  
one[,'issuing_divgov:opposing_next:next_term_divgov']=one$issuing_divgov*one$opposing_next*one$next_term_divgov
one[,'issuing_divgov:opposing_next:immutable']=one$issuing_divgov*one$opposing_next*one$immutable  
one[,'issuing_divgov:next_term_divgov:immutable']=one$issuing_divgov*one$next_term_divgov*one$immutable  
one[,'opposing_next:next_term_divgov:immutable']=one$opposing_next*one$next_term_divgov*one$immutable  
one[,'issuing_divgov:opposing_next:next_term_divgov:immutable']=one$issuing_divgov*one$opposing_next*one$next_term_divgov*one$immutable

# simulate predicted probabilities
zero=cbind(1,zero)
one=cbind(1,one)
p.zero=apply((1/(1+exp(-as.matrix(zero)%*% t(betas)))),1,as.vector) 
p_zero=apply(p.zero,1,mean)
p.one=apply((1/(1+exp(-as.matrix(one)%*% t(betas)))),1,as.vector) 
p_one=apply(p.one,1,mean)
diff=p_one-p_zero

diff_mutable_h3 <- diff

# results
h3_lb_a <- quantile(diff,.025)
h3_est_a <- mean(diff)
h3_ub_a <- quantile(diff,.975)

# effect of issuing overs under divided government for immutable orders

# negative case
zero=dt[complete.cases(dt[,vars[1:11]]),vars[1:11]]
zero[,'immutable']=1
zero[,'issuing_divgov']=0

zero[,'issuing_divgov:opposing_next']=zero$issuing_divgov*zero$opposing_next
zero[,'issuing_divgov:next_term_divgov']=zero$issuing_divgov*zero$next_term_divgov
zero[,'opposing_next:next_term_divgov']=zero$opposing_next*zero$next_term_divgov
zero[,'issuing_divgov:immutable']=zero$issuing_divgov*zero$immutable
zero[,'opposing_next:immutable']=zero$opposing_next*zero$immutable  
zero[,'next_term_divgov:immutable']=zero$next_term_divgov*zero$immutable  
zero[,'issuing_divgov:opposing_next:next_term_divgov']=zero$issuing_divgov*zero$opposing_next*zero$next_term_divgov
zero[,'issuing_divgov:opposing_next:immutable']=zero$issuing_divgov*zero$opposing_next*zero$immutable  
zero[,'issuing_divgov:next_term_divgov:immutable']=zero$issuing_divgov*zero$next_term_divgov*zero$immutable  
zero[,'opposing_next:next_term_divgov:immutable']=zero$opposing_next*zero$next_term_divgov*zero$immutable  
zero[,'issuing_divgov:opposing_next:next_term_divgov:immutable']=zero$issuing_divgov*zero$opposing_next*zero$next_term_divgov*zero$immutable

# positive case
one=dt[complete.cases(dt[,vars[1:11]]),vars[1:11]]
one[,'immutable']=1
one[,'issuing_divgov']=1

one[,'issuing_divgov:opposing_next']=one$issuing_divgov*one$opposing_next
one[,'issuing_divgov:next_term_divgov']=one$issuing_divgov*one$next_term_divgov
one[,'opposing_next:next_term_divgov']=one$opposing_next*one$next_term_divgov
one[,'issuing_divgov:immutable']=one$issuing_divgov*one$immutable
one[,'opposing_next:immutable']=one$opposing_next*one$immutable  
one[,'next_term_divgov:immutable']=one$next_term_divgov*one$immutable  
one[,'issuing_divgov:opposing_next:next_term_divgov']=one$issuing_divgov*one$opposing_next*one$next_term_divgov
one[,'issuing_divgov:opposing_next:immutable']=one$issuing_divgov*one$opposing_next*one$immutable  
one[,'issuing_divgov:next_term_divgov:immutable']=one$issuing_divgov*one$next_term_divgov*one$immutable  
one[,'opposing_next:next_term_divgov:immutable']=one$opposing_next*one$next_term_divgov*one$immutable  
one[,'issuing_divgov:opposing_next:next_term_divgov:immutable']=one$issuing_divgov*one$opposing_next*one$next_term_divgov*one$immutable

# simulate predicted probabilities
zero=cbind(1,zero)
one=cbind(1,one)
p.zero=apply((1/(1+exp(-as.matrix(zero)%*% t(betas)))),1,as.vector) 
p_zero=apply(p.zero,1,mean)
p.one=apply((1/(1+exp(-as.matrix(one)%*% t(betas)))),1,as.vector) 
p_one=apply(p.one,1,mean)
diff=p_one-p_zero

diff_immutable_h3 <- diff

# results
h3_lb_b <- quantile(diff,.025)
h3_est_b <- mean(diff)
h3_ub_b <- quantile(diff,.975)

t.test(diff_mutable_h3,diff_immutable_h3)

h3 <- diff_mutable_h3-diff_immutable_h3

# discussed in-text
quantile(h3,.025)
mean(h3)
quantile(h3,.975)

############
# FIGURE 1 #
############

# store results for main figure
results = rbind( c(h1_lb,h1_est,h1_ub,'(1) Immutability','All Orders'),
                 c(h2_lb_a,h2_est_a,h2_ub_a,'(2) Opposing Next','Mutable'),
                 c(h2_lb_b,h2_est_b,h2_ub_b,'(2) Opposing Next','Immutable'),
                 c(h3_lb_a,h3_est_a,h3_ub_a,'(3) Issued Div. Gov.','Mutable'),
                 c(h3_lb_b,h3_est_b,h3_ub_b,'(3) Issued Div. Gov.','Immutable'))

results=data.frame(results)

names(results)=c('lb.95','Marginal Change in Pr(Revocation Next Term)','ub.95','Expectations','Order Type')
for (i in 1:3) {results[,i]=as.numeric(results[,i])}

results$Expectations=factor(results$Expectations,levels=c('(3) Issued Div. Gov.','(2) Opposing Next','(1) Immutability'))

F1 <- ggplot(results) + 
  geom_hline(yintercept = 0,linetype='dashed')+
  geom_hline(yintercept = -.075,linetype='solid',color='gray95')+
  geom_hline(yintercept = -.05,linetype='solid',color='gray85')+
  geom_hline(yintercept = -.025,linetype='solid',color='gray75')+
  geom_hline(yintercept = .025,linetype='solid',color='gray75')+
  geom_hline(yintercept = .05,linetype='solid',color='gray85')+
  geom_hline(yintercept = .075,linetype='solid',color='gray95')+
  geom_point(aes(x = Expectations, y = `Marginal Change in Pr(Revocation Next Term)`, color=`Order Type`),size= 2,position = position_dodge(width = 1/2)) +
  geom_linerange(aes(x = Expectations, ymin = lb.95,ymax = ub.95,color=`Order Type`),lwd = 1,position = position_dodge2(width = 1/2)) +   
  coord_flip() + 
  scale_colour_manual(values = c('#000000','#75988d','#FF5C5C')) +
  theme(panel.background = element_blank(),axis.text.x = element_text(angle = 45))

############
# FIGURE 2 #
############

ts = summaryBy(immutable~year,data=dt,FUN=mean,keep.names=T)

names(ts)=c('Year','Immutable')
F2 <- ggplot(ts) + geom_point(aes(x=Year,y=Immutable)) +geom_smooth(aes(x=Year,y=Immutable),color='#FF5C5C')+
  scale_x_continuous(breaks=seq(1940,2020,15)) +
  theme_bw()



